Research on smoke simulation with vortex shedding

The Lagrangian vortex method has the advantage of producing highly detailed simulations of fluids such as turbulent smoke. However, this method has two problems: the construction of the velocity field from the vorticity field is inefficient, and handling the boundary condition is difficult. We present a pure Lagrangian vortex method, including a nested grid to accelerate the construction of the velocity field, and a novel boundary treatment method for the vorticity field. Based on a tree structure, the nested grid algorithm considerably improves the efficiency of the velocity computation while producing visual results that are comparable with the original flow. Based on the vortex-generating method, the least square method is used to compute the vorticity strength of the new vortex elements. Further, we consider the mutual influence between the generated vortex particles. We demonstrate our method’s benefits by using a vortex ring and various examples of interaction between the smoke and obstacles.


Introduction
Incompressible flows can be widely observed in our daily life, such as smoke from burning wood and the vortices generated by water flowing through obstacles. As a main research direction in fluid simulation, smoke simulation has drawn substantial attention in industry and academia. The effective use of limited computing resources to obtain a more realistic effect of smoke movement is still challenging.
Various categories of numerical methods were developed to explore different types of complex and changing smoke phenomena. In terms of computing variables, these methods are classified for solving velocity, vorticity, wave function, and so on. In terms of discretization forms, these methods are classified as Lagrangian methods, Eulerian methods, and Lagrangian-Eulerian hybrid methods. Because of their unique advantages, these methods can be used to numerically simulate various flows.
The Lagrangian vortex method has a unique advantage in the simulation of unsteady flows where the vortex structures play a leading role, such as the dynamical evolution of a forced shear layer [1], an impulsively started rigid body [2], and certain interacting coherent vortical structures [3,4]. First material elements themselves being discrete elements, the Lagrangian vortex method does not have to construct complicated Eulerian meshes for flows with complex geometries. Moreover, it significantly decreases numerical dissipations [5] by avoiding interpolating between grids with the physical quantities stored on material elements. Second, the Lagrangian vortex elements can be placed adaptively during the process of discretizing the vorticity field, with the number of vortex elements placed being proportional to the vorticity strength. In addition, the continuous incompressible velocity field can be built by using the Biot-Saviot(BS) law from the discretized vorticity field. Thus, a continuous solenoidal velocity field that fills the entire computational domain can be reconstructed with relatively fewer discretized vortex elements needed for dynamical evolution [6]. However, unlike the velocity-based Lagrangian method [7,8], the Lagrangian vortex method has not yet been applied widely in industry and academia for the following reasons. First, the interaction between vortex elements is non-local, leading to higher complexities of high-precision and high-speed algorithms for large-scale simulations [9]. Second, the boundary condition of the flow field is generally given based on the velocity field, making the conversion of the velocity to vorticity boundary conditions necessary [6,10]. This conversion could be nontrivial, particularly for moving boundaries.
To solve the first of the aforementioned problem, we use a nested grid to group the vortex particles based on their positions and treat the vortex particles in the same group as one super vortex particle in an approximate manner. When calculating the velocity at a given position, we traverse the super vortex particles instead of all the original vortex particles. Thus, we can reduce the number of traversals to decrease the computational cost. We use the vortex-generating method to handle the boundary condition for the second problem. Unlike the momentum form of the Navier-Stokes equation, the vorticity at the boundary cannot be directly set. The plane method [5,11] in classic airfoil theory is only suitable for simple non-slip boundary conditions. In this study, the generated vortex particles are placed on the boundary to satisfy the flow's non-through and non-slip boundary conditions. The vorticity of the generated vortex particles is implicitly determined by solving a linear system. Our contributions are summarized as follows: 1. A parallel nested grid method that groups the vortex particles to decrease the number of vortex particles traversed during the velocity computation, thus reducing the computational cost; [16], and vorticity confinement [17][18][19] approaches. A common feature of the Eulerian method is that the fluid domain should be discretized in advance. Hence, the Eulerian method is not suitable for smoke simulation in open spaces. Among the Lagrangian methods, the most widely studied methods are the smooth particle hydrodynamics (SPH) and Lagrangian vortex methods. The advantage of the SPH method is its high computational efficiency. This is because the SPH method does not have to solve the global Euler equation. However, its accuracy relies on the massive distribution of particles, and this is not suitable for smoke simulation in large spaces. To solve this problem, Ren et al. [20] further improved the performance of SPH smoke simulation by solving for the visible smoke without involving the transparent air.
The Lagrangian vortex method is the most efficient method for simulating complex flows with sparse date and minimal numerical dissipation. This method discretizes the vorticity field on the Lagrangian vortex element and implements a velocity field by using the BS law without the numerical dissipation that is associated with projection. The flow data can be stored on vortex particles [9,21,22], vortex filaments [3,[23][24][25][26], and vortex surfaces [27][28][29]. Based on the connectivity of the mesh, the vortex surfaces can be used to track the fluid surfaces. However, this discrete form cannot deal with the topological changes of vortex structures. The vortex filament discretizes the vorticity field on the closed curves, which naturally represents the vortex structure. A fundamental problem with the vortex filament method is the complicated curve repair operation to the topological changes in the filament, making large-scale parallel computing challenging. The vortex particle method, also known as the vortex blob method, generally replaces point vortices with certain vortex cores to remove the singularity in the kernel function. The vortex core has various shapes, such as sphere, ellipsoid, or small vortex sheet. Compared with the vortex filament and vortex surfaces, it has a lower computational overhead and better stability but with decreased accuracy.
One common problem of the Lagrangian vortex method is that the integration operation's time complexity, from the vorticity field to the velocity field, is O(N 2 ). To overcome this problem, some classic solutions were proposed, such as the Fast Multipole Method (FMM) [30], Particle-Particle Particle-Mesh (PPPM) [9], and Octree-Based Method (OBM) [22]. The FMM can reduce the integral operation's time complexity to O(N) and guarantee sufficient accuracy. However, this algorithm is considerably complex and difficult to deploy on new hardware, making it impossible to be widely used in industries. The PPPM method has a computational time complexity of O(NlogN) at the expense of partial accuracy and is relatively easy to implement. The OBM method uses a hierarchical tree data structure to group vortex particles, and the computational time complexity is O(NlogN). However, the tree structure that carries the vortex particles' spatial distribution information makes the parallelization of the algorithm challenging.
The boundary treatment is notoriously problematic in the Lagrangian vortex method. In computer graphics, researchers generally place vortex layers at the boundary to satisfy the non-slip and non-through boundary conditions. Based on this idea, Chorin [21] proposed a vortex-generating vortex method to deal with the boundary conditions of vortical flows. Park [31] obtained a reasonable vorticity distribution on the solid boundary by solving a system of 3K g equations, where K g is the number of generated vortex elements. Vines [32] reduced the size of the linear system to K g by applying only a non-through boundary condition to the normal of the solid, thereby solving the boundary integral problem. Appropriate boundary conditions can reflect the physical mechanism of vorticity generation on the solid. An SPH turbulence model can be used to deal with boundary conditions of vortex filaments [33].

Vortex dynamics
Vortex particle method. The Navier-Stokes equation can be written in two forms: the velocity (u)-pressure (p) form and the velocity (u)-vorticity (ω = r × u) form. The former focuses on tracking the change in momentum to solve the pressure, viscous force, or other external forces in the equation. In this study, we use the latter to simulate the smoke by using the following vorticity equation: where ν is the viscosity coefficient. Eq (1) expresses the incompressibility of the fluid. In Eq (2), the pressure term disappears. The left side of Eq (2) shows the change rate of vorticity under the Lagrangian viewpoint, which can be further expressed as @ ω/@t + (u � r)ω. Furthermore, (u � r)ω is the advection term. The right side of Eq (2) shows the stretching and diffusion terms, from left to right. The stretching term does not exist in 2D space. In 3D space, as a key element of turbulent motion, it can produce local intensification and reorientation of the vorticity. The diffusion term describes the fluid's vorticity diffusion due to friction, which can smooth the diversity of fluid movement and slow down the motion of flow. The diffusion term can be solved by employing the particle strength exchange (PSE) method [34], which is not discussed in this paper. Based on the Lagrangian viewpoint, the vorticity field is discretized on a set of vortex particles; this is consistent with the fact that the non-zero vorticity is generally concentrated in the flow trajectory [35]. The vorticity strength carried by the vortex particles can be computed by using following equation: where Γ i , ω i , and V i is the vorticity strength, vorticity, and volume of the i th vortex particle, respectively. Stretching. The methods used for solving the stretching term can be roughly classified into two categories. An explicit method is to discretize the velocity on the grid and compute the Jacobi matrix which is all partial derivatives of all components (x, y, z) of velocity u = (u, v, w). The accuracy of this method is considerably low because of the interpolation. Another approach is to capture the vortex structure's stretching by calculating the deformation of the geometric structure, such as the vortex ring or vortex sheet in the flow, which cannot be directly applied in the vortex particle method.
We use the vortex segment approach to solve the stretching term. To use vortex segments, which naturally represent the vortices' geometric structure, we build a bridge between the vortex particles and the vortex segments. Each vortex segment is a small cylinder of length h and constant circulation κ. The segment has the same direction as the vorticity strength and two endpoints, x l and x r . A vortex particle of vorticity strength Γ i can be translated to a vortex segment of unit direction We first convert the vortex particles into vortex segments and then update the ends of the vortex segments by applying the following equation: where u is the ends' velocity that can be calculated by using Eq (8).
Finally, we translate the vortex segments back into vortex particles with updated vorticity strengths. If we directly set the vortex particle with vorticity strengthΓ i nþ1 ¼ k i ðx nþ1 i;r À x nþ1 i;l Þ, even if the time step is small, the simulation diverges quickly. The reason is that, although the magnitude of velocity is small, the gradient of the velocity field can be large in a turbulent flow. This fact makes the velocity difference between the two ends significantly large, resulting in an unstable simulation. Therefore, we weightΓ i nþ1 and Γ nþ1 i to stabilize the simulation and obtain the final vorticity strength: where δ 2 (0, 1) is the interpolation coefficient. Notice that in our vortex segment approach, we smooth the update of the vorticity instead ofΓ i nþ1 to preserve sharp features as much as possible.
Velocity computation. In the Lagrangian vortex method, the velocity field can be obtained from the vorticity field by using the BS law: where p is the inquiring position, and d is the computational domain's dimension. In this study, the integration in Eq (6) is equivalent to the summation of all vortex particles 2D : uðpÞ ¼ 3D : uðpÞ ¼ where K v is the number of vortex particles, and x i and Γ i are the position and vorticity strength of the i th vortex particle, respectively. Moreover, Γ i is a scalar in 2D. The BS law has the following property: the formula's divergence-free velocity field is not unique. Consequently, the corresponding vorticity field is uniquely determined for a given velocity field. This property allows us to add a harmonic velocity field. However, the BS law has two problems. One is the expensive numerical integration. We introduce a solution to this problem in Section. Another problem is that vortex particle selfinduction can cause a singularity and numerical instability. Researchers have proposed certain methods to modify the BS law to circumvent this problem. Angelidis and Neyret replaced the BS kernel with a radial basis function [36], which is defined smoothly around the centroid. This makes it easier to analyze the integral.
Chorin [21] proposed the concept of a "vortex cluster," which is a small area with a high concentration of vorticity. The shape of the area can be circular, elliptical, or vortex sheet. The combination of the vortex clusters and different non-singularized kernel functions achieves different smoothing effects. In this study, an artificial smoothing parameter, σ, is added to Eqs (8) and (9) is our modified BS formula.
When the inquiring position, p, is far away from the vortex particles, K � 4π kp − x i k 3 . When p is close to the vortex particles, the properties expressed by Eq (9) are similar to the solid vortex core. Fig 1 shows the velocity profile induced by the modified BS law.

Fast summation
A major problem of the BS law is the time-consuming summation. As described by Eq (8), we need to traverse K v vortex particles when calculating the velocity of the inquiring position. This means that the time complexity of calculating the velocity of K v vortex particles is OðK 2 v Þ. When a phenomenon requires massive vortex particles to show the details, the computational cost is extremely high.
By using the BS law, we derive an approximate method for reducing the traverse number in summation. According to Eq (9), the velocity is inversely proportional to the square of the distance. This property shows that the distant vortex particles have less influence on the final result, whereas the closer vortex particles have a significant contribution. If the distance between the vortex particles is much shorter than their distance to the inquiring position, the distance between these vortex particles and the inquiring position can be regarded as approximately equal. By employing the approximation, we can obtain where K s is the number of vortex particles, p is the inquiring position, and � x is the weighted average position of the vortex particles.
By applying Eq (10), we find that a group of vortex particles have roughly the same contribution as a super vortex particle towards a distant inquiring position. Therefore, based on their positions, we can group the vortex particles and treat those in the same group as a super vortex particle. Thus, when calculating the velocity field, we only need to traverse the super vortex particles. This can reduce the computational overhead.
In this study, a nested grid algorithm that is similar to a tree structure is proposed to group vortex particles. The nested grid is a set of grid layers, G = {G i |i = 1, 2, . . ., K n }, with different resolutions and the same bounding box, where K n is the number of grid layers. The resolution relationship between grid layers G i and G i+1 is ðK iþ1 It means that each cell in grid layer G i corresponds exactly to the 2 d cells in G i+ 1 .
The consistency in space is convenient for program implementation. All grid cells are regarded as a super vortex particle, and its vorticity and position are calculated by using Eq (11).
where Γ s and x s are the vorticity strength and position of the super vortex particle; and K s , Γ i , and x i are the total number, vorticity strength, and position of the i th vortex particles in the cell, respectively. Based on the spatial consistency of different grid layers, we design the traversal manner of the super vortex particles as follows: 1. Starting from the first grid layer, G i=1 , traverse all the super vortex particles, C i¼1 s in G i = 1 , that is, the grid cells. Find all cells, C � s , that do not contain the inquiring position, p, and directly calculate their induced velocity. Then, find cell C � s that contains the inquiring position, p, and enter the next grid layer, G i=2 .
2. Starting from the second grid layer, G i=2 , traverse only the cell set, C i s , that corresponds to C � s in the upper grid layer. Find the cell set, C � s , that does not contain the inquiring position, p, in C i s and compute their induced velocities to accumulate the result with the upper grid layers. Then, find cell C � s containing the inquiring position, p, enter the next grid layer, G i+1 , and repeat the second step until the penultimate grid layer.
3. In the last grid layer, G K n , find the cell set corresponding to grid cell C � s in the upper grid layer and compute their induced velocity.
We use a 2D example to illustrate the operation of the nested grid. First, we group the vortex particles depending on their distribution in the computation domain, as shown in Fig 2(b). The vortex particles in the same group are regarded as a super vortex particle. Fig 2(c) shows the distribution of super vortex particles in different grid layers for calculating the velocity of p. Thus, the time complexity of calculating the velocity of p is reduced from O(K v ) to O(logK v ).
For games or real-time applications, performance and effects are more important than accuracy. Although our method loses some accuracy, it substantially improves computational efficiency and can obtain satisfactory visual results.

Boundary treatment
We use the vortex-generating method to deal with the boundary conditions for the Lagrangian vortex method. Owing to friction and viscosity, vortices are generated on the surface when the flow passes the obstacle. The vortex-generating method is consistent with this fact. The method can be divided into two steps: (1) compute the position and direction of the generated vortex particles; (2) calculate the magnitude of vorticity strength of the generated vortex particles. From simple to complex, we take the 2D static boundary example to introduce our method and then present the dynamic boundary example. Whether it is a static or dynamic boundary, we follow a principle when dealing with boundary conditions: the obstacle velocity equals the flow velocity.
First, we introduce the method of computing position x g and direction t g of the generated vortex particles. As shown in Fig 3, we sample some boundary points that are represented by the gray points with position x b . The number, K g , of boundary sampling points is the same as the number of generated vortex particles. To offset the slip and through velocity of the flow at the boundary sampling point, t g should be tangent to the velocity, u f , of the flow at the  boundary and the normal, n, of the solid surface. By normalizing the vorticity direction, we can obtain We place the generated vortex particle near the boundary along vector a ¼ a t g �u f kt g �u f k , as shown in Fig 3. As a result, x g = x b + a. The placement makes the generated vortex particle's induced velocity at the boundary opposite to u f . In our method, α is set artificially. The larger α is, the bigger the vorticity strength of the generated vortex particle, and vice versa.
Next, we use the least square method to compute the magnitude, Γ g , of the vorticity strength of the generated vortex particles. Assuming that the velocity of background flow is u 1 , we consider that the summation of the induced velocity, u g , of the generated vortex particles to the boundary and u f + u 1 is equal to the velocity, u b = 0, of the solid under ideal conditions. Therefore, the problem has been reduced to finding a minimizer of the total velocity.
We derive Eq (13) as which is the boundary velocity induced by the generated vortex particles with the unit vorticity strength vector. Eq (14) can be solved as Different from the static boundary, we need to compute the velocity, u b , of the solid for the dynamic boundary treatment. We also need to update the boundary sampling points in each time step, which can be computed by using Eq 16. The rigid body's motion is divided into two parts: the center translation and the rotation relative to the center. The position of the boundary sampling points is found as where x c ðtÞ and M c ðtÞ 2 R d�d are the center and rotation matrix of the rigid body, respectively. d is the spatial dimension, and ½ðx g ð0Þ À x c ð0Þ�M c ðtÞ T represents the rotation of the rigid body. Then, we take the time derivative of x b ðtÞ to determine the velocity, u b , of the rigid body at the boundary sampling point. Substituting Eq (16) and u b into Eq (15) yields the magnitude of the vorticity strength of the generated vortex particles at the dynamic boundary: where

Time integration scheme
Our time integration scheme can be summarized as follows. We first update the position of the vortex and tracer particles. In our method, tracer particles are used for visualization. Then, we handle the boundary condition and solve the stretching and diffusion terms. Finally, we construct the velocity field from the vorticity field by using the modified BS law and the background flow. In each time step, we update the states of the flow by using the following steps, as depicted in S1 Fig:   1. Advection. Calculate the induced velocity by using Eq (9) and the advect vortex and tracer particles by using the fourth-order Runge-Kutta method.

Parameter study
We conduct two sets of simulations with different parameter values to better understand the impact of the vortex particle's volume V and the stretching coefficient δ on the simulation. First, we simulate the rising smoke with three different values of V (as shown in Fig 4). In this example, the computational domain is approximately 20 � 20 � 20. We emit a vortex circle (consisting of a set of vortex particles) every 5-timesteps, whose normal towards the upper right. From left to right in Fig 4, V equal to 64υ 3 , 8υ 3 and υ 3 , respectively. Here, υ is 6.28 � 10 −3 .
To ensure the vorticity strength of the three flows are the same, the vortex particles' number of the vortex circle is 1000, 2000 and 4000 from left to right, respectively. We observe that the overall motion trends of the three cases are consistent, which shows that the volume of the vortex particle does not change the macroscopic behavior of the flow. In the fine scale, the difference of the vortices with different V are marked by the red circles. The vortices generated by the vortex particles with big V merge into larger vortices earlier. In contrast, the vortices in the flow with small V maintain longer. It is because the flow with small V has more vortex particles, which samples the vorticity field more uniformly. However, the increase in the number of vortex particles also brings more computational overhead. Therefore, we should balance the computational speed and animation quality according to the actual situation.
Second, we use three different values of δ in the example of colliding vortex rings (see Fig 5). We set a uniform grid inside the ink and place the vortex particles in the center of every grid. Except for the normal, the vorticity distribution of the two inks is the same, which is calculated by applying Eq (18). We observe that the vortex re-connection phenomenon can be obtained after the collision of rings in all cases, which is consistent with the result reported in [37]. It proves that the numerical damping caused by δ does not change the physical nature of the simulation. When δ = 1, the vorticity is not suppressed. The vortices after re-connection show irregular motion and instability. In the case with smaller δ, the flow after the collision is more stable. This proves that our stretching approach can stabilize the simulation. In addition, we compare our approach to a MacCormack solver [14], as shown in Fig 6. The simulation with a MacCormack solver can not reproduce such phenomena.

Nested grid results
Vortex ring. The simulation results of the octree grid [22], our nested grid and the direct summation are qualitatively compared by employing the vortex ring example. In the comparison, the initial conditions are set identically. The vorticity distribution of the vortex ring is approximated according to the empirical formula (18) [38]: where r � = r/R, R is the radius of the cross-section of the vortex ring, r is the distance from the vortex particle to the vortex axis, and K = 1/2exp(2)log(2). The vorticity is concentrated in the core of the vortex ring, and the vorticity direction is tangent to the vortex axis and perpendicular to the normal of the vortex ring. The nested grid is a five-layer grid with the configuration of ðK 5 r Þ 3 ¼ 64 3 . As depicted in Fig 7, our method's simulation results are very similar to those of the direct summation method. Although some details are lost on the fluid surface and the

Boundary treatment result
Kármán Vortex Street. Fig 9 shows the simulation of classic fluid phenomena, Kármán Vortex Street. In this example, we compare our method with the panel method [31]. At the boundary of the disc, we sample 200 panels with the panel method and 200 boundary points with our boundary treatment method. As shown in Fig 9, some vortex particles that enter the obstacle's interior are eliminated after advection. The generated vortex particles passing the obstacle produce a repeating pattern of swirling vortices, forming the Kármán vortex street phenomenon. Under the influence of viscosity, the vorticity strength of the vortex particles that are far away from the obstacle decreases slowly, and their velocities gradually become similar to the velocity of the background flow. However, the vortices in the panel method can not maintain for a long time, and finally merge together. This comparison showcase the effectiveness of our boundary treatment method. In addition, we simulate the vortex shedding behind two discs and an airfoil, as shown in Fig 10. We sample 200 boundary points around every disc and 400 boundary points around the airfoil. The airfoil's angle is 45˚. Compared with the disc, the flow passing around the airfoil presents a more irregular motion. The simulation shows that our method can handle the boundary of the obstacles with complex geometry.   In addition to the visual results, we provide the computation accuracy of the three cases presented in Fig 11. We define the ratio of velocity elimination as follows: where K g is the number of boundary sampling points, u f and u f+ are the velocity of fluid at the boundary sampling points before and after boundary treatment, respectively. u 1 is the velocity of the background flow, and u b is the velocity of the solid at the boundary sampling points. In the Kármán Vortex Street example, the solid is static and the velocity of the background flow is constant. Therefore, e can be simplified as e ¼ Fig 11 shows the ratio of the velocity elimination for the three cases. We can observe that our boundary treatment method can cancel out most of the slip and through velocities while dealing with regular obstacles. In the three cases, one disc achieves the highest accuracy, whose elimination ratio reaches approximately 83%. The accuracy of the airfoil is the lowest. We remark that our method is more suitable for dealing with the boundary of obstacles without sharp shapes.
Vortex shedding behind a static ball. Fig 12 indicates the realistic wakes behind a static ball in the rising smoke example. In this experiment, we emit a vortex circle every 10 time steps, consisting of 200 vortex particles. The vortex circle carries the tracer particles upward. When the static ball is passed, the vortex particles entering the ball are destroyed, and the collision point with the obstacle is the boundary sampling point. The radius of the ball is 0.5 unit lengths.  We construct the thin smoke to form the words 'PLOS One' and let it pass through the cylinders. The radius of the cylinder is 0.8 unit lengths. In this example, the boundary sampling point is computed in advance like the 2D Karman vortex street example. The vortex particles entering the cylinder are destroyed. The simulation results show that our method can produce dispersive thin smoke.
Vortex shedding behind moving objects. The first column in Fig 14 depicts the turbulent wake formed by the moving ball passing through the static smoke. Unlike in the 2D Karman vortex street's sampling manner, we place massive vortex particles regularly in the static smoke and treat the collision points between the vortex particles and the moving ball as the boundary sampling points. The vortex particles entering the obstacle are destroyed. This means that the number of vortex particles is constant in the simulation. The ball's radius is 0.5 unit lengths, and the ball moves at a speed of 2 unit-lengths/s. The second column in Fig 14 presents the turbulent wake formed by the moving and rotating bunny passing through the static smoke. In this example, the boundary sampling manner is the same as the moving ball. We use the signed distance field to calculate the bunny's normal and detect the collision between the vortex particles and the bunny. In addition to the translation speed, we introduce a rotation speed to the bunny to prove that our method can handle complex boundary conditions. Fig 14 shows that the bunny produces more details than the ball. The bunny's size is equivalent to that of the ball, the translation speed is 2 unit-lengths/s, and the rotation speed is 2 unit radians/s.

Performance
We parallelize most of our method's steps using OpenMP and implement a GPU version of the nested grid method using CUDA, which is the most computationally expensive part of our method. Our examples are performed on a laptop with a 6-core 2.20 GHz CPU and a NVIDIA GeForce GTX 1060 graphics card. Table 1 lists the timing statistics of various smoke examples.

Conclusions
We presented a pure Lagrangian vortex method to produce rich details for smoke animation. Our method is efficient and suitable for parallelism on a GPU, thereby indicating its potential to be applied in real-time applications. This framework has two essential components. A nested grid is used to group the vortex particles and effectively reduce the computational cost of the BS law. This is validated by using the vortex ring example. Further, we proposed a novel vortex-generating method that supports dynamic boundary conditions for handling complex obstacles. We simulated a series of complex gas phenomena, such as the Kármán Vortex Street, colliding vortex rings, and interaction of smoke and obstacles.
Future work: The nested grid loses some accuracy when computing the velocity. We can improve this by dividing a direct summation area in the finest grid layer, while the other layers remain unchanged. Moreover, the number of generated vortex particles determines the linear equation's size, making our method challenging in dealing with the interaction between smoke and large-scale obstacles. Reducing the linear equation's size or splitting the equation into multiple sub-equations could be beneficial. Finally, the two-way fluid rigid body interaction is a prominent research direction for future work.